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Abstract 

Experimental data regarding auxin and venation formation exist at both macroscopic and 
molecular scales, and we attempt to unify them into a comprehensive model for venation 
formation. We begin with a set of principles to guide an abstract model of venation formation, 
from which we show how patterns in plant development are related to the representation of 
global distance information locally as cellular-level signals. Venation formation, in particular, 
is a function of distances between cells and their locations. The first principle, that auxin 
is produced at a constant rate in all cells, leads to a (Poisson) reaction-diffusion equation. 
Equilibrium solutions uniquely codify information about distances, thereby providing cells with 
the signal to begin differentiation from ground to vascular. A uniform destruction hypothesis 
and scaling by cell size leads to a more biologically-relevant (Helmholtz) model, and simulations 
demonstrate its capability to predict leaf and root auxin distributions and venation patterns. 
The mathematical development is centered on properties of the distance map, and provides 
a mechanism by which global information about shape can be presented locally to individual 
cells. The principles provide the foundation for an elaboration of these models in a companion 
paper [13], and together they provide a framework for understanding organ- and plant-scale 
organization. 
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1 INTRODUCTION 



1 Introduction 

One of the principal tenants of biology is that no matter how large an organism becomes everything 
about it must ultimately have an explanation at the cellular level. Molecular biology goes even 
further by requiring an explanation on the level of chemical reactions. What chemical compounds 
and what reactions give rise to the intricate patterns of veins in leaves? Or to the pattern of 
specialized cells in the root of a plant? These are the types of questions that modern biology 
attempts to answer. And these same questions have prompted workers from the other sciences to 
join in. Physicists, mathematicians and computer scientists find the problems especially intriguing 
because of the need to explain how global patterns develop from local behaviors. Seen in this light, 
the problem becomes the search for a "plant geometry:" how cells determine where they are located 
with respect to other "special" cells, what those "special" cells are, and how cells behave once the 
information becomes available. 

Attempts to solve a version this problem, in which the notion of positional information is 
the focus, can be traced back to over a century ago [38, 39], but it was only within the last 
fifty years that mechanistic proposals were first submitted [41, 42, 18]. Of these, the idea of a 
diffusible morphogcn [37] has received the most attention because it captures complex measurable 
phenomena in a compact mathematical form. This so-called reaction-diffusion formulation requires 
specific knowledge of molecular interactions, but its typical form requires at least two chemical 
substances in order to explain a patterning phenomenon [20]. By contrast, we have shown [15] 
that a modification of the original formulation only requires one substance and already provides 
preliminary answers to the three main questions of plant geometry. 

This paper is the first of three, in which we develop these questions in further detail. The series 
looks for the simplest hypotheses that can explain patterning phenomena arising in vein forma- 
tion, facilitated transport of plant hormones, cell division and expansion, hormone concentration 
distributions, and others. We propose hypotheses about the local behavior of cells, such as how a 
hormone is produced, and analyze them mathematically to explain observed phenomena as well as 
to generate further hypotheses. Thus, some of our fundamental assumptions will be theoretically 
derived. The verification of their implications is presented through numerical simulations, which 
we demonstrate to afford unique interpretations. As a result, we develop a theory that explains 
how discrete systems - such as collections of cells - may compute a distance map in a variety of 
scenarios and represented by various interpretations of hormone concentrations. 
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2 CONSTANT PRODUCTION HYPOTHESIS: POISSON MODEL 



We begin, in this paper, by refining some of our assumptions in [15] about the biology of plants. 
To keep matters tractable, it is still necessary to include some abstraction for what otherwise might 
be considered signalling or other networks. We abstract "constant production", "proportional 
destruction", and c- vascular conversion "schema" in this paper. We provide the mathematical 
analysis of our earlier model, to which we refer as the Poisson Model, that proves our claims in 
[15]. Then, we extend the model by introducing a more biologically plausible assumption about the 
destruction of the signaling hormone auxin, which gives rise to the Helmholtz Model. Mathematical 
analysis of this formulation demonstrates that the properties of the Poisson model are kept and 
that it can explain even more experimental observations. 

In the next paper [13] we elaborate the schema into more biologically plausible mechanisms using 
different transport facilitators and the chemosmotic theory. We observe that it is remarkable that 
the Fickian transport and reaction diffusion equation developed here can provide this abstraction in 
such a manner that its main properties hold when more detailed facilitated transport is taken into 
account. Our goal throughout this series of papers [15, 12, 13, 14] is to formulate those abstract 
principles that can govern the qualitative properties exhibited at the systems level in plants. Such 
an approach is necessary, we believe, to organize into organ- and plant-scale syntheses the diversity 
of cellular and molecular mechanisms constantly being discovered. As we show through the series, 
an elaboration of the principles into increasing detail predicts non-linear and, at times, surprisingly 
delicate sequential developmental patterns. Without such a systems-level understanding one might 
be tempted to postulate a need for unnecessary genetic machinery. 



2 Constant Production Hypothesis: Poisson Model 

In [15], we proposed the Constant Production Hypothesis and argued that rich geometric informa- 
tion becomes available to cells in a single-substance reaction-diffusion model. 

Hypothesis 1 (CPH). Auxin is produced in all cells at the same constant rate. 

In this section, we recall the main consequence of this assumption and then prove it mathematically. 

2.1 Background 

A leaf is a collection of cells. We distinguish between ground cells, those that give rise to all others, 
and vascular cells, those that comprise the venation pattern. We focus on early leaf development 



Draft of May 27, 2009 [15:20] 



2.2 Poisson Model 



2 CONSTANT PRODUCTION HYPOTHESIS: POISSON MODEL 



Poisson Model 
Definitions 

Ground Cell: Diffusion coefficient Dg. 

C- Vascular Cell: At least one interface has diff. coef. By > Dg. 

Cell Functions (Program) 

CFl: Produce substance s at constant rate K. 

CF2: Measure c and Ac through interfaces. 

CF3: Diffuse s through interfaces. 

CF4: When Ac > r through interface /, change its diffusion coefficient to Dy. 

Tabic 1: Poisson model definitions. The mechanism for changing the diffusion coefRcent in CF4 
will be elaborated in [13]. 

and concentrate on signals sufficient to initiate the cascade of events that change ground cells into 
vascular cells within an expanding arcolc. To keep matters tractable, a cell will be referred to 
as c-vascular (cascade vascular) immediately after this cascade is initiated. The role of this c- 
vascular abstraction is to summarize the increasingly elaborate cascade of genetic expression and 
transcription regulation that is being uncovered; see [31]. The sub-collection of c-vascular cells 
may be thought of as an early pre-pattern from which veins derive. Ground cells have (essentially) 
homogeneous characteristics and areoles are delimited by more developed c-vascular (or mature 
vascular) cells. Instead of assuming the pre-pattern is predefined, our model establishes how it 
emerges from local operations. We refer to both membranes and cell walls together as cell interfaces 
and assume that they act as a single membrane. 

2.2 Poisson Model 

Each cell performs the basic functions listed in Table 1 independently and simultaneously. Under 
these assumptions, then, cell functions CFl and CF3 determine the equation governing the distri- 
bution of s in the areole. They define how the substance is produced and transported for both 

ground and c-vascular cells. The latter evacuate the hormone much faster so we assume that the 
boundary of the areole may be thought of as a sink for s, i.e. it is essentially kept at a constant 
level. Therefore, the temporal change of the concentration inside a region depends on how much 
diffuses in or out of a cell plus how much is created; in symbols. 



Ct = DV^c+K, 



(1) 
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2.2 Poisson Model 2 CONSTANT PRODUCTION HYPOTHESIS: POISSON MODEL 

where D is the diffusion constant of ground cells, V^c = Cxx + Cyy is the Laplacian of concentration 
over cell position, and K is as in CFl. This is a reaction-diffusion equation which has a steady- 
state: after a sufficiently long time, the dynamical system is well approximated by the c such 
that Ct = (see Fig. 1). Observe that those cells which are further from the boundary have 
higher concentrations. In fact, the concentration profile is qualitatively similar to that of the 
function assigning to each cell the shortest distance to a (c-)vascular cell — the so-called distance 
transform [3]. 

When Ct = 0, Eq. 1 becomes a standard Poisson equation. Given our boundary conditions 
(c = at veins), there is a unique c satisfying it [16]. Prom this, we calculate: 

Result 1. Consider an areole and suppose that P is a ground cell which is furthest from the c- 
vascular boundary. Let Q he a c-vascular cell which is closest to P and denote by L the distance 
between P and Q. Then 

(a) c(P) is proportional to ^L"^; 

(b) the change in c at the interface of Q nearest to P is proportional to j^L. 

(c) Ac is largest at an interface of the c-vascular boundary, larger than for any ground cell, and 
is proportional to ^L. 

Therefore, using part (a) and CP2, a cell may determine if it has become further than L units 
from the closest c-vascular (supply) cell by measuring its concentration. More must be done, 
however, to guarantee that the developing vascular network is connected, and utilizing the difference 
in concentration accomplishes this. 

Result 1(b) asserts that Ac at the venation is directly proportional to L/D and does not depend 
on the value of c. Moreover, it also gives the direction toward the furthest cell. This is sufficient to 
show that mechanisms for new strand creation should adhere to the following schema: 

Schema 1. Let Dj be the diffusion constant across an interface I and Ac be the concentration 
difference through I. Then increase Dj to a higher value when Ac > a^Lo. (a is a constant of 
proportionality.) Alternatively, the flux cf) = DjAc = uKLq may be employed. 

An illustration of this mechanism is shown in Fig. 1. 
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Left Middle Right 



Figure 1: Hormone concentration inside an areole. (left) A rectangular domain (artificial areole) 
is illustrated with a boundary of c- vascular cells. Assuming c- vascular cells are much more efficient 
at transporting s, the boundary may be taken as a sink and c is governed by Eq. 1. (middle) c at 
near steady-state, Cf « 0. Also shown are the values of c and Ac along a path (in black) across the 
areole. Notice how the concentration peaks at cells furthest from the veins, while Ac peaks near 
the vein, (right) (C) Concentration, (D) magnitude of gradient, (E) gradient vector field. Observe 
how the gradient vectors point toward largest concentration increase. 

2.3 Analysis of the Poisson Model 
2.3.1 Definitions and Background: Geometry 

A collection of ground cells surrounded by c-vascular cells is called an areole. Our technical result 
will assume that an areole is a discretization of a continuous portion of which we call a shape. 

Definition 2. A shape is any subset G which is the closure of a bounded open set and has a 
boundary dfl consisting of finitely many smooth curves. 

A point Q G dCl is concave if for any line £ locally tangent to Q there is an open ball B^{Q) such 
that B,{Q)r]lr\{n-dn) = B,{Q)r]l-{Q} (i.e. the line segment is inside n). If Se(Q)n^nf^ C dn, 
then Q is a corwea; point. The boundary has concave (convex) curvature^ at concave (convex) points. 

Definition 3. Let f2 be a shape and P G M^. The Euclidean distance function on f2, denoted £n^ 
is 

^We allow infinite curvature. 
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The boundary support of P, denoted bsupp(P; 9f2), is 

bsupp(P; dQ) = {Q e 9fi : ||P - Q|| = £n{P)} ■ 

The medial axis of fl, denoted MA(f2), is the set of points P which have two or more closest 
points on the boundary, i.e. 

MA{n) ^{PGn-. |bsupp(P; dn)\ > 2} . 

where |bsupp(P; (?ri)| denotes the cardinality of the set. Note that MA(il) does not have to be 
restricted to the shape fi and is well-defined on all of M^. Hence, there is an interior medial axis 
and an exterior one. Here, we will only be concerned with the interior one. 

Theorem 4. Let fl be a shape and P G fl. Suppose |bsupp(P; 9f2)| = 1 and pick the unique 
Q e bsupp(P;9f}). Then, 

(a) P ^ dCl implies 

where Q — P is the vector from Q G dfl to P. 

(b) Ifdfl is C'' at Q, then VSq is C'' at P. 

Proof. Part (a) is due to [17, 4.8(3)] and (b) is a consequence of the more general result by [23] (see 
also [25]). □ 

Corollary 5. £q(P) is smooth at P Gfl- MA(f2). 

Proof. Immediate from Theorem 4(b) and the definition of shape. □ 
Theorem 6. Let Cl be a shape. Then 

(i) MA(f2) has no interior, i.e. it is thin. 

(ii) MA(r2) consists of a finite number of connected piece-wise smooth curves. 

(Hi) if P G MA(f2), Q e bsupp(P;9f2) and C is the center of curvature for dfl at Q, then 
\\P - Q\\ < \\C - Q\\ whenever \\C - Pjj < jjC - Qjj. 
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Figure 2: Examples of distance map and medial axis (see [33]). top: negative distance map —£n, 
CENTER: level sets of £o, bottom: medial axis computed as in [11]. 
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Proof. Part (i) is shown in [26] and in [5]; (ii) is treated in detail by [7]. Part (iii) asserts that if a 
medial axis point is inside the circle of curvature of a point in its boundary support, then it cannot 
be further than the center of curvature. □ 

Theorem 7. Let Q be a shape. There is a unique confl such that c = ondQ and V^c = —K/D. 

Proof. See [8, p. 246] or Theorem 4.3 of [19] for a more general statement and proof. Also see 
[21]. □ 

Theorem 8 (Divergence). Let dBs{P) be a circle or radius s centered at P G M."^ , J\f the inner 
normals. Then 

V'^c{P) = lim / (Vc, A^) ds. 
Proof See p. 151 in [40]. □ 
Definition 9. The 0-notation for asymptotic behavior of a function is defined as: 

Q{g{n)) = {/(n) : 3ci,C2,no positive s.t. Vn > no,0 < c-ig{n) < f{n) < C2/(n)} . 

2.3.2 Statement of Result 

Result 1 is based on the following theorem. 

Theorem 10. Let Cl be a shape and c : — > M the unique function satisfying c{x,y) = on 
{x,y) e dCl and 

V^c=-^ . (2) 

Suppose P is such that £n{P) = L = sup^fo and Q € bsupp(P; d^). Suppose the smallest 
concave curvature radius is pL with p> 0. Then, 

(a) c(0) G e(L2), 

(b) i,L < |Vc| < f L^, 

(c) SUP90 |Vc| > supo_an |Vc| 

If an areole is regarded as a discretization of a shape fl, then the discrete approximation behaves 
as stated in Theorem 10. Observe that Vc at dfl is perpendicular to the boundary because c = 
there; hence, Vc{Q) points in the direction of P according to Theorem 4. 
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2.3.3 Organization of the Proof 

Parts (a) and (b) of Theorem 10 follow from Lemma 20. The idea of the proof is to find appropriate 
bounding functions, one from below v and another u from above, that sandwich the unique solution c 
and that take the same values at the boundary. Thus, v < c< u everywhere and \Vv\ < | Vc| < | Vu| 
on points where v = c = u, i.e. the boundary. Since the value of c must be the same at the boundary, 
its gradient there must be perpendicular to the boundary which gives the direction as claimed in 
Result 1(b). Lemma 14 and Lemma 18 give the lower bound and upper bound constructions and 
Lemma 20 collects them. 

Part (c) of Theorem 10 is necessary to show that the c-vascular strand creation process is 
well defined. Result 1(c) is the non-technical version of this claim which is stated more precisely 
in Lemma 21. The proof is based on the idea that the boundary may be seen as evolving by 
considering level sets of c, i.e. points 7c(, where c{x, y) = cq. The gradient must be perpendicular to 
this level set and the solution of Eq. 2 inside it follows the same constraints as the shapes on which 
the problem is defined. We may move the level set curve so that the point on which is on the 
gradient curve initiated at the point Q of maximal gradient on 70 touches Q for small enough cq. 
Knowing that the solution on the smaller shape must be strictly smaller than on the original shape 
shows that the maximum gradient magnitude must be strictly decreasing as the curve evolves. This 
is true for all curves, including the evolved ones (i.e. 7^^ for ci > Cq), so the gradient in the interior 
of the shape must be lower than the maximum on the boundary. 

2.3.4 The Proof 

We begin with Lemma 11 which will be used (indirectly) in most of the proofs that follow. It states 
that a discretization of the dynamic process will always move the concentration values in the same 
direction (up or down) if this direction is locally the same for all discrete points. This fact will be 

used to prove the next result. Lemma 12, which states that the equilibrium solution over a shape 
completely contained in another shape will be bounded above by the solution over the larger shape. 
This holds even if the initializing function is not smooth. 

Lemma 11. Let Ct = DV^c+ K he approximated on a square lattice by pi and its four neighbors 




where h is the lattice spacing. Suppose that Ct < (>)0 



everywhere on the domain of definition at time to. Then 
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(a) c+ TCt will also satisfy the inequality ifO<T< and 

(b) the discrete dynamics with such r make c decrease (increase) monotonically everywhere. 

Proof. Part b) follows directly from a). Let a = Djl? and Ap = X]j=i 4 '^('^j) ~4c(p), so Ct(fo,f>) = 
aAp + K. After the time step r the approximation to becomes: 

Ct(io + T) = £'V2(c+TCt(to))+i^ 

= ^ ^X^(c(%) + rct(to,nj))-4(c(pi)+TCt(io,Pj))^ + 
= a ["aj,, +r ^^ct(io,nj) - 4ct(io,j)i) j j + -ff 

= a I^Ap. + TQ A„ , - 4Ap^ j j + if 

4 

= aAp. (1 - 4rQ!) + ara ^ A„^. + 

Now, all of the ah.n^-\-K < and aAp. < by assumption. Also, < 4ra < 1. So choosing the 
largest A„j and replacing for the other three bounds the value above (since A„^ < 0), shows that r 
is used in a linear interpolation between two non-positive numbers. This finishes the claim. □ 

Lemma 12. Let Ut = DV'^u + K = over fl with u = on dfl. Iffl'cfl and DV'^c +K = Q on 
Q! with c = Q on dCl' then c< u onCl'. Iffl^fl', then c <u everywhere on fl' — dfl'. 

Proof. If 16 = over fl, then ut = K > and, by Lemma lib, u > in the interior fl — dfl after 
any non-zero time step. Thus, at equilibrium, u will satisfy the dynamics everywhere on fl' except 
possibly on dfl' where it may need to be lower because it the boundary conditions. If a boundary 
point is lowered in the discretization of the problem, any neighbor will see its V^u decrease and 
strictly decrease if the neighbor is not moved. This holds at any time step and will affect all points 
after sufficiently long time because they are all connected. If fl fl' , then u is not a solution (u > 
somewhere on dfl') and the dynamics on fl' will strictly monotonically lower it everywhere. □ 

Now we show what the solution looks like in one dimension. Lemma 13, and then we turn to 
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two special shapes: the circle (Lemma 14) and the open doughnut (Lemma 15). These shapes will 
be instrumental in providing the lower and upper bounds needed later on. 

Lemma 13. Suppose the domain is the segment [0, L], ci£)(0) = and Cid'{L) = 0. Then the 
solution to Eq. 2 is 

Proof. By inspection since the solution is unique: d"^ /{dr)^[ciD{r)] = —K/D. □ 

Lemma 14 (Disc). Let the shape be a circle of radius L centered at the origin. Suppose c{x,y) = 
on the boundary where a;^ + = . Then the solution to Eq. 2 is 



and 



\\/c{x,y) \ = -^^Va;2 + 2/2 



Proof. By inspection since the solution is unique. We see that c^x = Cyy = ~\K/D, so Eq. 2 is 
satisfied. □ 

Lemma 15 (Open Doughnut). Let < I < I + L be the radii of two circles centered at the origin. 
Suppose c{x,y) = on the boundary x^ + y'^ = P and Vc{x,y) = for x^ + y^ = {I + L)^. Then 
the solution to Eq. 2 is 

c{x,y) = CD (^Vx"^ + 2/^) 

where 

cd{t) = CD{r; UL) = ^ Q - r') + kn (^) {I + Lf^ (3) 
Further, c{x, y)>0 for l'^ < x'^ + y^ < {I + Lf. 

Proof. By inspection since the solution is unique. Since c = at the inner boundary and Vc 
points radially toward the outer boundary, the values of c are increasing radially in P < x"^ + y"^ < 
{l + Lf. □ 

The next two results (Lemma 16 and Lemma 17) are technical assertions used in the proof of 

the Upper Bound Lemma (Lemma 18). This is the last result needed to prove Lemma 20 and, 
therefore, parts (a) and (b) of Theorem 10. 
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Lemma 16. Let Cl be a shape P G fl, and Q G bsupp(P; dfl). Let S be the circle of curvature of 
fl at Q. Then £q{P') = £^{P') + 0{e^) for \\P — P'\\ = e, and over a circular region R of radius 
e centered at P 

lim ^ / (^CD{Sn)^£n,M) ds = lim ^ [ ( -^cd{Ss)^S,:,^) ds 
e^o 7r£2 Jgji \ dr I e^o wE^ J \ dr I 

Proof. Write the second order approximation = f s + 0[e^^ and Vfn = Vfs + In a 

circular neighborhood i?, the limit becomes; 

Now, ^Cd{£t. + Oie^)) = ^cd{£t) + Oie") by inspection of ^C£)(r), which gives: 



_ 1 9 



(|:Cb(£:e)V£:5: + |:Co(£:e)0(£2) + 0(£3)V£:e + Die')) 



1 d 

lim — ^cu(£'i;)V£'e 



□ 



Lemma 17. Suppose the shape Q is i/ie disc as in Lemma 14 with radius l + L and c is the solution. 
Then, 

u{x,y) = cd{1 + L - \/x'^ + 2/2) 
satisfies V^m < V^c = —K/D for all points except the center. 

Proof. Write f = u — c where c is the solution for the disc from Lemma 14. Letting R{x, y) = 
\/x'^ +y'^ 

f{R{x,y)) = cd{1 + L-R)-\c^d{1 + L-R) 

= %{^i + \{lr.{^){l + Lf-{l + L-R){l + L) 

A direct calculation shows that 

2 IK l + L 
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which demonstrates that /(r) < everywhere and linir^o fir) = — oo in the center of the disc. 
Therefore, u = V ^ + V^c < -K/D because V^c = -K/D. □ 

Lemma 18 (Upper Bound). Let the conditions of Theorem 10 hold. Define u{x,y) = 2cd{1 + 
£n{x, y)) with I = pL. Then u> c. 

Proof. We shall show that any discretization u with spacing h < ho (for some ho > 0) of u will 
satisfy DV^u + K < 0. Thus, Lemma 11 shows that a dynamical process initialized with u will 

decrease u everywhere with each time step and, by Theorem 7, it should converge to c. Hence, 
c< u. 

First, we treat non-medial axis points. Let P = {x,y) £ Q (not on the medial axis) and Q e dQ 
which is closest to P, i.e. \\P — = S{P). Suppose dO. near Q is approximated by the circle of 
curvature at Q. Thus, Lemma 16 applies and, by the Divergence Theorem 8, V^u(P) is the same as 
if the boundary were a circle at Q. Following Theorem 6 MA2, there are three cases: (a) P outside 
the circle, (b) P inside, and (c) the circle has infinite radius - it is a line segment. Lemma 15 shows 
that V^w(P) = —2K/D < —K/D which takes care of (a), and (b) is covered by Lemma 17. If the 
boundary is locally a straight line, then 

= ^^cn{l + r) = -^{l + {l + Lfy) (4) 

with I < r < I + L and r is in the direction of the gradient. So, V^u(P) < —K/D for all 
Per2-MA(f2). 

Now suppose that the of u on f2 are sampled on a discrete square lattice with spacing /i > 0. 
There, V^u(p) is approximated by the formula for Ap in the proof of Lemma 11. The error of the 
approximation is 0{h^) (see [1]). Thus, from the above, h may be chosen so that Ap < —K/D for 
P gQ. further than h from MA(f]). 

Let to be the distance function from the circle of curvature at the boundary point corre- 
sponding to P ^ MA(n). Hence, if ||P' - P|| = /i, then (P') =£q{P')+sp^ where Ep^ = 0{h^). 
Set 

Ep = sup |ep/| and s= sup ep' 
\\P'-P\\=h PGn-MA(n) 

where < e < ft, for small enough h. So, choose such an h and define 



Uh{x,y) = 2cd{1 + £n{x,y);l,L + 2h) 
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and notice that we may refine the grid (i.e. choose h smaller) and the above properties will still 
hold. Thus, refine h if necessary to make < —K/D on shape points further than h from the 
medial axis. Refine it further to A^ < —K/D on an open doughnut with I inner radius and L + 2h 
outer radius. Make sure that h is small enough so that c_D(r + h) — 2c£)(r) + CD{r ~ h) < —K/D 
(this is needed in the tangent line construction below), which is possible because of Eq. 4. 

Now we show that this also makes A^ < —K/D for points on the medial axis and those closer 
than h from it. Pick such a P and let consider Q £ bsupp(P; dil). If Q is concave, then approximate 
the boundary by its circle of curvature and look at = 2cd{1 + £i:{x, y); l,L + 2h) . A neighbor 
N of P used in A^ satisfies £j:{N) > £q{N) (because Q is concave and the difference is no more 
than s. Hence, u^{N) > Uh{N) since cd is increasing until I + L + 2h. Further, £i:{P) = £n{P) 
because the circle of curvature touches dCl at Q. Hence, Ap** < Ap^ < —K/D because P is not a 
medial axis point for E. 

If, on the other hand, there is a concave Q G bsupp(P; dCl), then instead of the circle of curvature 
we may take the tangent line £p at Q define S^p exactly similarly to S-^^ above. Refine h so that 
any point P' for which ||P — P'|| = h is closest to a point Q' £ £p that lies outside the fl or 
on af2.2. Thus, as before, £ip{N) > £n{N) for any N near P, i.e. - P|| = h. Therefore, 
A^" < Ap^ < -K/D. 

Finally, Lemma 12 shows that Uh > c from which we conclude that u> c since lim/i^o Uh — > 

u. □ 

Remark 19. The function u{x, y) need not be smooth on f2. In fact, it will fail to have fist derivatives 
on the medial axis of most shapes. 

Lemma 20. Let the conditions of Theorem 10 hold. Then c = Q{L'^). Further, if P is the center 
of the largest inscribed circle and Q a point on the boundary of the shape and the circle, then Vc{Q) 
points toward P and 

—L < Vc < —L ^ . 

2D - ^ ' - D p 

Proof. The Disc Lemma 14 gives the lower bound function and the Upper Bound Lemma 18 the 
rest. The gradient points in the direction of the normal to the boundary because c = on dfl. 
The magnitude follows from a simple calculation of d/dr[cD{r)] a,t r = pL (see the Open Doughnut 
Lemma for the definition of Coir)). □ 
^This must be possible since Q is convex. 
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2.3 Analysis of the PoissonmdamSTANT PRODUCTION HYPOTHESIS: POISSON MODEL 

Finally, the following result completes the proof of Theorem 10. 

Lemma 21 (Decreasing Gradient). Let c satisfy the Poisson equation (Eq. 2) on Q. and c = Q on 
dn. Then 

M = sup |Vc| = sup |Vc| 

and 

\Vc{x,y)\<M, {x,y)Gn-dn. 

Proof. Let jco = {{^^ y) & ^ ■ c{x, y) = cq}. A number < cq < supfj c must exist since c > on 
O — dn by Lemma lib. Let the shape il' be defined by G such that c{x,y) > cq. Hence, 

fi' C O and 0^0'. The boundary d^l' = -yco is regular, so there is a unique v satisfying Eq. 2 on 
ft' with u = on dQ' . Thus, v = c — cq and Vv = Vc on Q'. 

Thus, 7o = dCl and is connected for small enough cq (because is the closure of an open 
set). Further, if cq < £o for some sq > 0, then is a smooth curve because V-y^^c = on 7cq, 
c is at least twice differentiable, and Vc ^ when taken over Cl on points of 70. In fact, Vc is 
perpendicular to the curve 70 = dfl. Let Q € 70 be such that Vc{Q) = supg^ |Vc|. Let /3 be the 
integral curve segment starting at /3(0) = Q with tangents in the direction of Vc and such that 
/3(1) € jca- Since Vc is perpendicular to 70, /3(1) will be the closest point to 70 from 7cj, for small 
enough cq. Thus, may be translated so that /3(1) touches Q ensuring that jco is completely 
contained in Q; denote this translated curve by j'^^. 

The solution v' to Eq. 2 on 7^^ and its interior must be the translated v. By Lemma 12 < c 
everywhere except on 70 r\j'ca (e.g. at Q) where v' = c. Hence, |Vt;'(Q)| < |V'u((5)|. Therefore, 
|Vu(Q)| is strictly decreasing in the direction of Vu(Q). 




(a) (b) (c) 

Figure 4: Level sets of c. (a,b) Rectangular shape, (c) Areole from [31]. 

□ 
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3 PROPORTIONAL DESTRUCTION HYPOTHESIS: HELMHOLTZ MODEL 



3 Proportional Destruction Hypothesis: Helmholtz Model 

3.1 Background 

In [15] we noted that new veins may emerge simultaneously in areoles of drastically different sizes. 
We argued that the most parsimonious explanation of this phenomenon is that auxin is produced 
at a rate that is constant for each cell, regardless of that cell's size. But we did not develop the 
question of auxin destruction beyond assuming that veins drain the hormone in such a way that it 
is effectively depleted from the areole. A more natural assumption emerges from considering recent 
molecular work. 

Auxin is involved in a variety of processes that take place in plant cells at all times, so a 
portion of the free lAA is constantly being depleted. For example, it has recently been shown 
[22, 10] that the early response genes are activated by auxin through an increased degradation of 
promoter inhibitors. Thus, auxin binds to a TIRl protein which is instrumental in tagging inhibitor 
proteins (Aux/IAA) for degradation: the larger the concentration of auxin, the more effective the 
degradation. In addition, the hormone "is readily conjugated to a wide variety of larger molecules, 
rendering it inactive. Indeed, the majority of lAA in the plant is in the form of inactive conjugates. 
Auxin conjugation and catabolism can therefore decrease active auxin levels." [35, p. 853]. For 
these reasons, we propose the Proportional Destruction Hypothesis: 

Hypothesis 2 (Proportional Destruction). Free auxin levels are constantly being depleted at a rate 
proportional to the available hormone levels. 

The constant of proportionality will be denoted by a and referred to as the destruction constant. 

3.2 Helmholtz Model Definition 

The Proportional Destruction Hypothesis leads to the updated formulation of our model shown in 
Table 2. We shall assume that the only mode of auxin transport is diffusion according to Pick's Law. 
Our constant production of K (kg-s~^) is replaced by a speed of increase in concentration equal to 
K/S{i) where S{i) is the size of cell i. The proportional destruction of auxin implies a decrease 
in concentration given by —am{i)/S{i) = —ac{i) with c{i) denoting concentration. Therefore, the 
manner in which the auxin concentration changes with time within each cell can be written as: 

gi=DV^c+--ac. (5) 
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3.3 Analysis of 8heBR[M10RTMmAL DESTRUCTION HYPOTHESIS: HELMHOLTZ MODEL 



Helmholtz Model 
Definitions 

Ground Cell: Diffusion coefRcient Dg-, production p = K/S. 

C- Vascular Cell: At least one interface has diff. coef. D^> Dg\ production p = K/S. 
Cell Functions (Program) 

CFl: Produce substance s at constant rate K and destroy it at a rate ac. 
CF2: Measure c and Ac through interfaces. 
CF3: Diffuse s through interfaces. 

CF4: When Ac > r through interface /, change its diffusion coefficient to Dy. 

Table 2: Helmholtz model. Terms definitions: K is the pcr-ccU substance production rate 
(mass/time); S is the size of the cell (volume); p is the per- volume production rate; are diffu- 
sion coefficients related to the permeability of cell interfaces; c is the concentration of the substance 
inside the cell, a is the destruction constant, and Ac(7) is the difference of concentration through 
the interface /; r is a threshold. 

The equilibrium of this dynamical system, when ct = 0, is therefore described by an inhomoge- 
neous Helmholtz equation, which is why we refer to this formulation as the Helmholtz Model. 

3.3 Analysis of the Helmholtz Model 

It turns out that this relatively small modification to the Poisson model greatly improves the model's 
descriptive power. We now establish the mathematical results that make this claim concrete before 
we turn to the discussion of biological experiments in the next section. We demonstrate that 

the same type of distance information is captured by the Helmholtz model in Proposition 22 and 
Proposition 24, but that it can be obtained in more than one fashion if the destruction constant is 
small. If this constant is large, we show that a different kind of distance becomes available to the 
discrete system - one given by the logarithm of concentration. 

Proposition 22. Consider the dynamical system 



Suppose that it acts over a domain fl which a shape as in Theorem, 10 and on which we impose 
a zero-flux boundary condition (Neumann). Let p^i : Q R. Then the following holds. 

(a) If a > 0, then lim^^^ c= c^ for a unique steady-state c„. 

(h) Let a = and R = J Pndfi/ Jdfi he the average production. Then limj^oo Q = R and c 
converges to Ca + est. whenever R = 0. Further, Vcq is unique even when R^ 0. 



do 
dt 



= DV^c + pq, - ac. 



(6) 
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3.3 Analysis of 8heBR[M10RTMmAL DESTRUCTION HYPOTHESIS: HELMHOLTZ MODEL 



(c) If A, B G R, then the transformation pn i— > ApQ + aB induces a unique transformation of the 
steady state Ca Aca + B and vice versa. It follows that the gradient of Ca is only affected 



Remark 23. In part (c), if the destruction term is not linear, e.g. ac + (3c^, then the gradient might 
be affected by B as well. 

Proof. Parts (a) and (b). To show existence we prove that the dynamical system achieves Cj = 0. 
Consider the dynamical system Cu = D'V^Ct —acf The boundary conditions arc inherited: since no 
flux goes through the boundary, there must be no change of concentration in time, i.e. Vcj • n = 
on d^. The unique solution of this system is Ct = 0. 

To prove uniqueness, suppose Ui and U2 both satisfy the equation given the boundary conditions 
and Ct = 0. Thus 



which gives rise to DV'^v — av = where v = ui — U2 and Vv • n = where n is the normal to 
the boundary. Since v is elliptic and a > 0, v vanishes everywhere and uniqueness follows (see [8, 
p. 329 and 321]). The same reference shows that if a = 0, then this uniqueness is up to an additive 
constant u = Ui + est; that is, only Vu is unique. 

Now to show the convergence in (b) whenever R = 0, note that Cu = DV^Ct assuming a = 0. 
This has a steady-state s.t. q — est. everywhere. Also, J Ci = J pjjdO which shows that C( = R. 

Part (c). Let Ca satisfy Eq. 6 for Cf = and a production function Pq'^ . Then, Z?V^c„ — aCa = 
—Pq^. Suppose c = Aca + B satisfles the equation for some pQ. Since this c is unique, the following 
verification proves the claim. 



We can relate the steady-state solution of Eq. 6 with small a to the steady-state solution of 



if Ay^ 1: Vca 1-^ 



AVCa. 



DV^ui + pn - aui = £>V^U2 + Pn - au^ 



DSI'^c - ac= -pn 

DW^{Ac^ + B)- a{Ac^ + B) = -pn 
ADV'^Ca - Aaca - aB = -pn 
A{DS/^Ca - aCa) = -pn + aB 
M-P^n^) = -pa + aB 
pn = A{p^^^) + aB 



The other direction is derived similarly and the result follows. 



□ 
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the dynamical system in the previous section. In fact, we now show that there are conditions 
under which the two systems are similar even though the boundary conditions are different. The 
key difference is that before we assumed a constant value for c at the boundary whereas now we 
only assume no flow through the boundary. Fig. 5 illustrates the correspondence in 1-D, and the 
following proposition makes the claim in 2-D precise. 




10 20 30 40 50 

Figure 1 . One dimensional simulations. Gradient of concentration at steady-state 

(a) (b) 

Figure 5: Illustration of the steady-state in a 1-D system with two domains of production values 
and small a. (a) Uniform diffusion coefficients, (b) Smaller diffusion coefficients on the left, top: 
Concentration at equilibrium (blue curve) and production curve (green curve) on y-axis against 
cell position on x-axis. BOTTOM: Gradient of concentration (or difference in concentration Ac). 
Notice that Ac attains a maximum where the production (equivalently, cell size) changes abruptly. 
Observe also that the location of this maximum has similar properties to the location of a sink in 
1-D: e.g., the portion in (a) where K ^ A = 1 has a parabolic concentration profile with a maximum 
where the Ac is lowest and a minimum where the Ac is largest. 



Proposition 24. Let be a shape with two components = ^IqUCIi such that floHfli = dfla. Let 
Dq and Di be the diffusion coefficients inside f^o and fii respectively. If J^^ Podw + /^^ Podv — 
and po(rio) — K J^^ dv > 0, then 

lim Ca — ck 

Do/Di^O 

where ck satisfies Theorem 10 for the shape CIq by setting c/<(9rio) = 0. 

Proof. The convergence of the system derives from Proposition 22(b). As Dq/Di — > the relative 
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speed of diffusion in fii increases to infinity. Thus, the concentration over fii will tend to a constant 
and, consequently, so will c{dflo) = c{flor\fli). The conditions of Theorem 10 are therefore satisfied 
and the claim follows. □ 

Note, however, that the destruction constant must be sufficiently small in order to obtain a good 
correspondence. But the value of a has a more important role. A strictly positive a implies that 
there is a maximal distance beyond which information cannot travel. Assuming that measurements 
that cells can perform have a limited precision, our next result shows that the contribution of cell 
A to the concentration at cell B will be negligible whenever A is sufficiently far from B. The larger 
the value of a, the shorter this distance needs to be. 

Theorem 25. Consider the dynamical system in Eq. 6 and assume the conditions on the domain 
as in Proposition 22. Let G{r;a) = exp(— r^/(2a-^)) and consider the Caussian convolution kernel 
Gi = (aV27r)G(v/a;2 + ?/2;c7), a = 1/a. Then 

lim (Gi * pn) = Ca 

where -k denotes convolution and c^ is as in Proposition 22. 

Proof. By inspection, c* = Gi ★ pa satisfies Eq. 6 as a — > oo. 

Suppose that we have a convolution kernel G2 for which cr is a function of a and such that 
aG2 = 1 and (j{a) — > as a — > 00. Therefore, G2* pn ^ S so S — a{G2*S) — > 0. We now show 

that Gi has this property, and that V^(G2 ★ S*) — » which shows that 

lim £»V^c* + pa-ac* =0 

a— ^00 

and proves the claim. 

The extrema of Gxx = dG^ / dx^ are at r = — 0, — (t\/3, (T^/i. The values are Gxx{^) = 

-1/ct2G(0) = -1/0-2, and G^:,{a^/2,) = 2/0-2 exp(-3/2). Hence, choosing Gi = a^j^f^G implies 
that sup|aGf/aa;2| = 0{a) and that ^Gi = 0{a) (because /^^ G^^x^ + a) = 1/(0-^277). 
Thus, setting a = 1/a, we have that V'^{Gi * pn) = ct — > and that / aGi = 1. The claim now 
follows. □ 

Corollary 26. Suppose a shape has two components = iloUili such that ilof^^i = dilo- If 
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Pni^o) = and pn(f2i) = 1, then 



a- 



lim 



oo 



logCa(Q) - logC, 



-a 



= est. > 



where Q € dClo and Sq^ is the distance function on CIq. 

3.4 Experimental Support of the Helmholtz Model 

The principal biological support of our Poisson model is an indirect one: the model produces 
patterns that are similar to patterns observed in nature. But it does not, for example, predict 
the concentrations of auxin in any measurable fashion. By contrast, our Helmholtz model does 
make such predictions and some experimental data is available. Ljung et al. [24, p. 466 and Fig. 1] 
"observed an inverse correlation between leaf size and lAA concentration that was independent of 
growth conditions and developmental stage." They measured the proportion of hormone mass to 
total leaf mass piaa (with units pg-mg^^ ) in leaves of different weight, W , and obtained data that 
can be described well by a function piaa — AW^^ where A is a constant and x ranges between 
0.72 and 0.98. This suggests that all cells may be producing auxin at the same constant rate if the 
hormone is depleted proportionally to its concentration. We reason as follows. 

Let m{i) denote the mass of auxin in cell i and consider how this quantity changes with time as 
auxin is produced, destroyed and transported to and from neighboring cells. The production is a 
constant, K, and the destruction, as we argued above, is proportional to the available amounts so 
this rate of change can be expressed as mt{i) = K — am{i) + Transport. Assuming that the leaf is 
detached, as it is prior to measurement, the transport term only moves the hormone inside the leaf 
but does not contribute to either a total increase or a total decrease. Therefore, the rate of change of 
auxin mass in the whole leaf is Mt = ^mt{i) = ^ K — m{i). The equilibrium of this system, 
when Mt = 0, describes well the state of the leaf during measurement because the leaf is small (most 
leaves are less than 10 mg) and the time needed to perform the manipulations — during which these 
dynamics apply exactly — is therefore sufficiently long to shift the distribution of the attached leaf 
to this equilibrium. Consequently, the constant production hypothesis predicts a total auxin mass 
of around Mjaa = nK/a for a leaf with n cells. The quantity reported in the Uterature, however, is 
an auxin-to-leaf weight ratio which we calculate to be piaa = Mjaa/W. Ljung et al. [24] plot such 
ratios for four groups of leaves of different sizes against leaf weight. Group members are selected 
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Figure 6: Inverse correlation between leaf size and lAA concentration. After Figure 1 of [24]. SD: 
short day; LD: long day. Original caption follows. lAA levels in Arabidopsis leaves. The lAA 
concentration was measured in (a) leaves 8-20 from six plants grown for six weeks under SD (b) 
leaves 1-12 from six plants grown for 4.5 weeks under SD (c) leaves 1-7 from three plants grown 
for 3 weeks under SD and (d) leaves 1-8 from 10 plants grown for 16 days under LD. The data are 
presented as loglog"^ plots of the lAA concentration in individual leaves versus leaf weight. 

according to the order of leaf emergence (phyllotaxis) and are organized as follows: the first group 
contains samples from leaf numbers 8-20, the second from leaves 1-12, the third from leaves 1-7, 
and the fourth from leaves 1-8. Leaves within each of these ranges have a fairly similar final shape 
and size (Ref. [36]) from which we infer that all leaves in the same group have roughly the same 
final number of cells. Therefore, since most data points are obtained after cell division has ceased 
and the cell numbers have stabilized, our analysis suggests that Mjaa is roughly the same for all 
samples in the same group and that only W differs. Theoretically, then, we expect the curves to be 
described by a function piaa = AW^^ with A a constant and a; = 1 which is in good agreement 
with the experimentally derived values of x w 0.72, 0.74, 0.80, 0.98, Fig. 6. 

3.5 Further Predictions of the Helmholtz Model 
3.5.1 Vein Formation 

The new (Helmholtz) formulation of the model preserves the distance information available locally 
to cells under appropriate conditions. Thus, if a is small and there are consistent differences in cell 
size between ground cells and c-vascular cells, then the distribution of auxin in an areole encodes 
size information just as it did in the Poisson case (details in Proposition 24). For example, if ground 
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Figure 7: Our schema for the elaboration of new c-vascular strands is unaffected by the relative 
sizes of boundary cells and interior cells — a consistent difference is enough. However, when all cells 
perform the functions CFl and CF3 from Table 2 the equilibrium concentration of s is substantially 
different for the two types of configuration. (A— C) C-vascular cells are larger than ground cells; 
(E— G) depiction of c-vascular cells (green squares) as new strands form; (H— J) c-vascular cells 
smaller than ground cells. (D,K) Arrows show direction of auxin flow. Here a = 0.01 in all 
simulations. 

cells are smaller than c-vascular cells (as in Fig. 7 A), then the same qualitative distribution of auxin 
is obtained as in [15]. The program in Table 2 then creates new strands as before. Such relative 
cell sizes are observed in the early emergence of c-vascular networks (see p. 21 in [28]), but in more 
mature tissues it is the ground cells that are larger (e.g. the procambium in Fig. 2 of [27] or the 
mature vein cells compared to others in [30, p. 234] or [29, p. 460]). Our model accommodates this 
second possibility as well. If boundary cells are smaller than interior cells (as in Fig. 7T), then the 
hormone distribution will be inverted — higher concentration on the boundary than in the interior — 
but the differences in concentration between neighboring cells will follow the same qualitative rules 
as in the first relative size configuration, albeit with an opposite sign (details in Proposition 22). 
Therefore, the same program can produce new vascular strands and the strands that it produces 
will be exactly the same in both configurations. 

Fig. 7 shows a simulation of new strands created according to Table 2 in a hypothetical areole 
for each of the two combinations of relative cell sizes. Notice that the two new c-vascular strands 
connect in the middle of the areole. They do so by draining s in opposite directions so that the 
cell where the two strands eventually meet does not have a well defined polarity — it is effectively 
bipolar. Thus, our model predicts multi-polar cells whenever loops of veins form. Recently at least 
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the bipolar case has been reported for Arabidopsis [32, Fig. 2], Fig. 8 compares our predictions to 
the empirical observations. 

However, our theory explains this phenomenon only under the assumption that vascular cells 
are larger than ground cells. If these relative sizes are reversed, then we should expect bipolar 
cells to form but the definition needs to be revised. In effect, the relevant cells would not have 
carriers facilitating export; instead, neighboring cells from opposite sides of a bipolar cell both 
exhibit facilitated transport toward the bipolar cell. Our theory predicts that such configurations 
would arise if new vascular strands are created in more mature organs, extending an existing vein 
as opposed to existing procambium (e.g. a tertiary vein stemming from a primary may be a good 
candidate). 

3.5.2 Auxin Distribution in Arabidopsis Roots 

Roots have a much simpler cell size distribution than leaves, and reliable estimates are easier to 
obtain. There even exist detailed 3-D models of root tips of Arabidopsis, but we shall only use 
2-D slices to test our theory in those organs. This type of data is representative of the full 3-D 
organ because roots are radially symmetric: rotating the 2-D slice about its long axis yields a 
good approximation the complete root. Along the same axis, three regions are distinguished: (1) a 
division zone (DZ) near the tip, (2) an elongation zone immediately adjacent to the DZ, and (3) a 
maturation or differentiation zone ([29, p. 436], [30]). On average, cells are smallest in the division 
zone, followed by slightly larger ones in the elongation zone, and then by the largest cells in the 
maturation zone. 

This schematic configuration is depicted in Fig. 9^ and forms the setup of the first type of root 
simulations. Each cell is represented by a rectangular box containing a single dot inside, the unit of 
auxin production independent of cell size. Note that the hormone distribution obtained predicted in 
this fashion (Fig. 9B,E) qualitatively agrees with the empirical measurements (Fig. 9C,D) reported 
by Bhalerao et al. [2]. Their data come from slices of untreated plants and consist of the average 
concentration of lAA as a function of distance from the root tip. The technique for measuring 
auxin levels gives good concentration estimates but has poor spatial resolution, so our comparison 
is restricted to the overall shape of the curve. The use of staining techniques, on the other hand, 
promises higher resolution but only provides qualitative information. Thus, we suspect that there 
is a peak of auxin concentration near the root tip, as shown in Fig. IOC, but do not know its height. 
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This observation may be sufficiently explained by the geometry of the organ, as our first schematic 
simulation suggests, so we turn to a real root specimen to test this claim. 

Our next simulation uses a manually traced 2-D slice [6, Fig. lA] to predict the shape of auxin 
distribution inside that root. The setup is redrawn in Fig. 10^4. Notice that the simulation result 

(in part B) reproduces that peak. This is a robust feature of the model as the only external cue 
has been cell size: the diffusion coefficients are all equal, the per-cell auxin production is constant. 

4 Numerical Simulations 
4.1 Background 

Here we provide background definitions and results separated in two major categories: geometry 
and graph theory. The first category is needed to prove the results about the steady-state of certain 
dynamical processes, and the second category is used to prove that a discretization of these processes 
is well behaved. 

A graph G = {V, E) is a combinatorial structure that consists of three things: a set of vertices 
(or nodes) V, a set of edges E, and an incidence relation. The edges describe which vertices are 
connected and the incidence relation attributes an order to each connection. For example, the edge 
e = («,j) says that node i is connected to node j and that it "starts" at i and "ends" at j. The 
number of edges of which vertex j is an end is called the degree of i. Each edge may take a value, 
called a weight, and these values can be recorded in an adjacency matrix A of size n x n where 
n = \V\. So, Aij is the value for the edge Note that the matrix will be symmetric if the 

direction of the edges does not matter, i.e. Aij = Aji. This is the undirected case, which will 
be discussed in this paper and describes the (Helmholtz) dynamics of auxin concentration. In the 
following paper [13], however, directed graphs will be needed and Aij = —Aji in some cases. 

A matrix can be described by its eigenvectors and eigenvalues. A vector v is called an eigenvector 
(or characteristic vector) of the matrix A iff Av — Xv for some number A, which is the corresponding 
eigenvalue (or characteristic value). 

The graph Laplacian L{G) of a graph G is given by L{G) = D{G)—A{G) where D is the diagonal 
degree matrix {Da =degree of i) and A is the adjacency matrix. The following is a standard result 
(e.g. see [4]). 
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Theorem 27. Let G be a connected graph on n vertices and denote by Xo < Xi < ■ ■■ < A„_i the 
eigenvecotrs of L = L{G), the graph Laplacian. Then: 

(a) All Xi are real; 

(b) Ao = 0, with eigenvector 1 = [1, 1, • • • , 1]-^/ 

(c) Ai > 0; and 

(d) the eigenvectors of L span M". 

A matrix A is called positive definite if v^Av > (where is the transpose of v) for all non-zero 
vectors v, and it is semi-definite if some non-zero vector w exists such that v^Av = 0. Note that 
L{G) is positive semi-definite. 

The determinant of an n x n matrix A is given by the following formula: 

n 

det{A)= sgn{a)l[A,^„(^i) 

where a is a permutation on n elements and sgn{a) is the sign of the permutation: positive if the 
permutation can be produced by an even number of element exchanges and negative otherwise. A 
matrix A is invertible - i.e. A~^ exists such that A~^A = 7rf - iff det{A) ^ 0. The determinant is 
equal to the product of eigenvalues, det{A) = Xi, so L{G) is not invertible. 

4.2 Discrete Simulations 

In this section we prove that appropriate discretizations of the continuous equations may be solved 
numerically. Specifically, we show that Eq. 1 and Eq. 5 converge to a unique solution as t — > oo 
and show how to obtain this solution without iteration. We shall not assume anything about the 
dimensionality of the space which contains the cells, only that they are connected. 

Suppose there are n cells in a conglomerate Cl where each cell shares an interface with at least 
one other cell and that the conglomerate is connected in this fashion. Labeling each cell with 
a number from 1 to n, suppose that each cell i contains a hormone at concentration c{i). The 
interfaces allow this hormone to diffuse following Pick's law, so assuming the diffusion constant 
through an interface between cell i and cell j is Dij , then the flow into cell i through each existing 
interface is given by Dij{c{j) — c{i)). Thus, the difi^usion of the hormone through the conglomerate 



Draft of May 27, 2009 [15:20] 



29 



4.2 Discrete Simulations 



4 NUMERICAL SIMULATIONS 



may be described by a matrix —L applied to the vector c of concentrations, where 




Notice that L is the graph Laplacian for the graph where each cell is a node and the edge weights 
are the diffusion coefficients Dij. 

Now suppose that the hormone concentration in cell i is somehow maintained at a fixed level 
c{i) = cf, refer to such a cell as a sink. Then the diffusion process in or out of cell i has no effect on 
its concentration c{i), but neighboring cells c{j) will be affected by c{i). Hence, the new matrix M^- 
describing the diffusion process looks exactly like —L except for the rows corresponding to sinks; 
if i is a sink then Ma = 1 and Mjj = for j ^ i. This means that det{M) = det{M"^), and that 
M is no longer symmetric, because any neighbor j of i which is not a sink will induce Mji = Dji. 
However, the sub-matrix M"* of M with rows and columns corresponding to all cells which are not 
sinks is symmetric. In fact, the next lemma follows immediately. 

Lemma 28. Let be a conglomerate ofn — Ug + Ur cells of which Us are sinks. Label the sinks 1 to 
Ug and the rest with + 1 to n. Let M"" — M^n^ + 1 : n, + 1 : n) be the x n^. sub-matrix and 
—L be the graph Laplacian of the sub-graph G^"*^ of non-sink cells. Let Cij by an Ur x diagonal 
matrix such that for each cell i in neighboring sinks sj we have Cu = '^g. Dig^ . Then, for a 

suitable labeling, 



If we also assume that the hormone is destroyed in each non-sink cell i at a rate proportional 
to the concentration in cell i, then the sub- matrix becomes 



where / is the n^. x identity matrix. Observe that M is invertible if and only if M"* is invertible 
because det{M) = det{M'^^). 

Lemma 29. // either a > Q or a > Q and there is at least one sink cell, then the matrix M is 
invertible and has strictly negative eigenvalues. 

Proof. As observed above, it suffices to show that det{M"^) ^ 0. We shall use the following result 



= -L-C . 



= -L-C-aI 



(7) 
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and apply it to -M""* from Eq. 7. 

Lemma 30. Let D = {dij} be an n x n diagonal matrix with da > with at least one da > 0. 
Let L he the graph Laplacian of a connected graph on n vertices. Then M = L + D is symmetric 
positive definite. It follows that M is invertible and has strictly positive eigenvalues. 

Proof. Prom Theorem 27 we see that L is positive semi-definite. So, if v e ffi" and v 7^ 0, then 
v^Lv > with equality reached only for v = si where s 7^ 0. Similarly, D is positive semi-definite 

because Dv = davf and all terms are non-negative. Further, (sl)^ D{sl) = dus'^ > 
since > and at least one entry in D is strictly positive. Finally, if v G M" and v 7^ 0, then 

v'^Mv = v^(L + D)v = v^Lv + v^L>v > . 



This finishes the proof of Lemma 30. 



□ 



Now, let 2? = C -I- al. If there is at least one sink, then at least one Cu > 0; if a > then all 
Da > 0. Hence Lemma 30 applies and completes the proof of Lemma 29. □ 

To complete the discretization of the continuous process, suppose that pn{i) = K/S{i), with 
units {mass * time~^ * volume~^), is the rate at which cell i produces the hormone. Setting S{i) 
to be the size of cell i, we obtain the following discrete dynamics: 



(8) 



dc/dt 



where At is the time step. We wish to show that, given a small enough At, c will converge to a 
unique value. 

Observe that the update rule in Eq. 8 may represented as 0^*+^^*^ = i/c^*^ by writing 



/ c(*+^*)(l) ^ 



n{t+At) 



(2) 



c(*+^*)(n) 
1 

c(t+At) 



(At)M + 1 



Aipn(l) 

Atpn(2) 

1 



4t) 



(1) 



J V 



cW{2) 

c(*)(n) 
1 



£(*) 
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where / is the n x n identity matrix. The convergence of the process now reduces to showing that 
jjk^io) converges as fc — > oo. Notice that U has the block form of the following three matrices 



A Va 




B Vb 




C V, 












1 




1 




1 



where A, B and C are n x n matrices; v^, V(, and Vc are n x 1 vectors; and is a 1 x n vector. 
The product of two such matrices preserves the block form; e.g. Uc = UaUb by setting C = AB 
and Vc = Avb + Va- Therefore, by induction on k, the blocks of Uc = U^ must be C = and 

Theorem 31. Suppose that each cell i has size S{i), produces a hormone at a rate pn{i) (ind 
destroys the hormone at a rate ac{i) with a>0. If a> or there is at least one sink cell, then for 
a sufficiently small At the discrete process in Eq. 8 will converge to c* = —M~^{p^/S). 

Proof. The conditions of Lemma 29 apply, so M has strictly negative eigenvalues and exists. 
Choose < At < 1/A where A is the largest (in absolute value) eigenvalue of M. Thus A = 
AtM + 1 will have eigenvalues < |Ai| < 1; it follows that limfe=(x) A'' = 0. Now, from the above, 
Vc = Eto ^'v„ = {A'' - I){A - = {A'' - 7)(AiM)-iva and the claim follows. □ 

This result demonstrates that an iterative process will indeed converge — assuming perfect arith- 
metic operations — but it also shows that the equilibrium can be computed much more efficiently. 
It suffices to solve the linear system c* = —M~^{pq/S). The system is well behaved numerically 
whenever a is sufficiently large, because the condition number of this matrix is roughly equal to 
the largest degree of the graph, times D divided by a; see Dahlquist and Bjorck [9] for a discussion 
of matrix condition numbers. 

4.3 Geometric Domain Definition 

In this section we outline how the domains — representing leaves, roots, etc. — are defined geometri- 
cally and then converted into the graph representation discussed in the previous section. Ultimately, 
the geometry of the domains should correspond to and be comparable to the geometry of real plant 
tissues. Thus, we define the domains by manipulating images of those tissues. Both the cell size 
and the cell neighbors (i.e. the topology of the graph) are computed from an image. 
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In this paper we adopted a pixel-based approach whereby the organ is drawn as a digital image 
and the color of each pixel encodes some information: whether the pixel is part of the domain or 
not, the value of the production function p, whether the pixel is a sink or part of the vein pattern, 
etc. The natural connectivity of pixels on a square grid — four or eight neighbors — then defines the 
topology of the graph. The final diffusion matrix is built after defining the diffusion constants for 
each pair of pixel colors. 

But this representation also allows us to define a cell by using multiple pixels. Fig. 9 shows an 
example in which a cell consists of several black pixels — representing the interior — surrounded by 
green pixels — representing the cell walls. None produce auxin except for a single dark-green pixel 
in the middle of the interior pixels. Thus, each cell produces auxin at the same rate (the rate of the 
center pixel) but cells may have different sizes. Moreover, the diffusion coefficient in the interior of 
the cell may be different from the diffusion coefficient through the cell wall. Our usual assumption 
is that the interior diffusion coefficient is much larger. 



5 Conclusion 

We have developed the foundations for a theory of how global information about shape is related 
to the distance transform, and how several of the essential properties of this distance transform 
can be computed by a simple reaction-diffusion equation. The model has its roots in our earlier 
Constant Production Hypothesis, and is based on a computational abstraction that all cells behave 
according to the same rules. Most importantly, it provides a mechanism that illustrates how "hot 
spots" of concentration can develop from structural conditions rather than differential production 
induced by an explicit developmental program. 

The explicit assumption about hormone depletion — the Proportional Destruction Hypothesis — 
greatly increased the scope of our earlier model [15]. We showed that there are at least two 
additional ways in which a distance map becomes locally available to a group of cells, and that 
testable predictions ensue. And although the available data is insufficient to compare the predictions 
to the actual numbers, the qualitative trend is accurately reproduced and explicit measurements 
are suggested by updated model. The analysis and assumptions of Section 3.4, in effect, describe 
further experiments to test the theory. 

The simulations in this paper, which involved detailed anatomical considerations, show the 
power of such calculations. That an auxin concentration peak emerged properly near the root tip 
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illustrates their role is sufRcency rather than necessity. 

Nevertheless, our model is still too abstract to be deemed biological. In particular, it is well 
known that diffusion is not the only transport mechanism responsible for auxin flow. Active, or 
at least facilitated, transport carriers are known to exist, which our current formulation does not 
consider. That is a topic of our companion paper. For now, wc remark that the Fickian transport 
and reaction diffusion equation developed here can provide an abstraction in such a manner that 
its main properties hold when more detailed facilitated transport is taken into account. 
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Figure 8: Illustration of the constant production hypothesis and cartoon model for vein formation. 
(A) Consider a square areole as in Fig. 7A-D in which c-vascular cells are larger than the interior 
ground cells. Auxin diffuses faster between c-vascular cells than any other type. We show the 
equilibrium concentration distribution. Note that it is minimal nearest the veins and maximal at the 
center; i.e., it varies with the distance to the nearest vein. Arrows along a one-dimensional cut (black 
line through the center) illustrate the flow of auxin along this line from high- to low-concentration 
pixels. (B) Concentration along the cut in A illustrating maximum at center. Magnitude of 
gradient (concentration difference between cells) varies "inversely" and peaks at the veins with 
a value proportional to vein distance to the center. This suggests a schema : differentiate from 
ground to vascular when gradient magnitude is large (equivalently: when central cells are far 
from veins). Once a cell begins to differentiate, it clears auxin more efficiently thereby causing 
adjacent cells to differentiate, until new veins are formed (see D below). (C) PINl expression 
in Arabidopsis. Red star (*) denotes a bipolar cell. (D) Cartoon mechanism of vein formation 
suggested by Scarpella et al. [32] based on measurements as in C. Note that this realizes precisely 
the schema in A. (E) Illustration of the veins formed in an areole according to the schema in A and 
developed in [15]. Note the bipolar flow at single cells predicted by the model and reported in [32]; 
compare with C. Colors denote magnitude of gradient and arrows show the direction of vascular 
strand formation, opposite to flow. (F) Cell outlines for the areole modeled in E. Note the large 
vascular cells. Figure credits: A, B, i? from [15]; C, D from Fig. 2 and Fig. 7, resp. of [31]; Ffrom 
[28]. 
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Figure 9: Simulations in a schematic root. 



E 




(A) Setup inspired by known regions of 
plant roots [29, p. 434-6]. Cells are rectan- 
gular boxes which increase in size from the 
tip toward the shoot. A single dot inside 
each cell represents one unit of hormone 
production; green boundaries (cell inter- 
faces) have the same diffusion coefficients 
everywhere. (B) Result of simulation. 
Concentration at steady-state through the 
horizontal red line in A. (C,D) Measure- 
ments of auxin concentration in sampled 
root tissues (either 5mm or 1mm cylinders) 
reported by Bhalerao et al. [2, Figs. 2 and 
4]. (E) Result of simulation: concentra- 
tion profile over the whole domain. 
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Figure 10: Root simulations from a manually traced Arabidopsis specimen from Swarup et al. [34, 
Fig. 1]. (A) Traced root used for the simulation. Cell size was determined by computing the 
area; neighbor relations and interface area were obtained by computing the length of a shared 
boundary. Diffusion coefficients and per-cell hormone production are the same for all cells. (B,C) 
Result of simulation. Note the predicted peak of concentration (arrows). (D) Stained root from 
Casimiro et al. [6, Fig.lA] showing a peak at the tip (arrow). Our model predicts this peak even 
though all parameters are kept uniform. 
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